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Abstract 



A dynamical mean-field approximation (DMA) previously proposed by the 
present author [H. Hasegawa, Phys. Rev E 67, 041903 (2003)] has been 
extended to ensembles described by a general noisy spiking neuron model. 
Ensembles of A^-unit neurons, each of which is expressed by coupled K- 
dimensional differential equations (DEs), are assumed to be subject to spa- 
tially correlated white noises. The original i^A^-dimensional stochastic DEs 
have been replaced by K{K + 2)-dimensional deterministic DEs expressed 
in terms of means and the second-order moments of local and global vari- 
ables: the fourth-order contributions are taken into account by the Gaussian 
decoupling approximation. Our DMA has been applied to an ensemble of 
Hodgkin-Huxley (HH) neurons {K = 4), for which effects of the noise, the 
coupling strength and the ensemble size on the response to a single-spike in- 
put have been investigated. Results calculated by DMA theory are in good 
agreement with those obtained by direct simulations. 
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I. INTRODUCTION 



It is well known that a small cluster of cortex may contain thousands of similar neurons. 
Each neuron which receives spikes from hundreds of other neurons, generates spikes prop- 
agating along the axon towards synapses exciting neurons in the next stage. Dynamics of 
an individual neuron with voltagc-dcpcndcnt ionic channels can be described by Hodgkin- 
Huxley-type (HH) model [1], or by reduced, simplified neuron models such as integrate- and- 
fire (IF), FitzHugh-Nagumo (FN) [2,3] and Hindmarsh-Rose (HR) models [4]. Although the 
response of a single neuron in vitro is rather accurate, that in vivo is not reliable [5]. This 
is due to noisy environment in living brains, where various kinds of noises are reported to 
be ubiquitous (for a review sec [6]). In recent years, the popiilation of neuron ensembles 
has been recognized to play important roles in the information transmission (pooling effect) 
[7]- [12]. Then it is necessary for us to theoretically investigate high-dimensional, stochastic 
differential equations (DEs) describing the large-scale noisy neuron ensemble. In order to 
make our discussion concrete, let us consider ensembles consisting of A'"-unit neurons, each of 
which is described by X-dimensional coupled DEs: for example, K —1, 2, 3 and 4 for IF, FN, 
HR and HH neuron models, respectively. Dynamics of such neuron ensembles, expressed 
by ii'A'"-dimensional stochastic DEs, has been so far investigated with the use of the two 
approaches: (i) direct simulations and (ii) analytical methods such as Fokker-Planck equa- 
tion (FPE) and the moment method. Simulations have been made for large-scale networks 
mostly consisting of IF neurons. Since the CPU time to simulate networks by conventional 
methods is proportional to N^, it is rather difficult to simulate realistic neuron clusters in 
spite of recent computer development. In FPE dynamics of neuron ensembles is described 
by the population activity. Although FPE is the powerful method formally applicable to 
the case of arbitrary K and A^" [13], actual calculations have been made mostly for A" = oo 
ensembles of a A' = 1 model with the use of the mean- field and/or diffusion approximations 
[14]. Similar population density approaches have been recently developed for a large-scale 
neuronal clusters [15,16]. The moment method initiated by Rodriguez and Tuckwell (RT) 
has been applied to single FN [17,18] and HH neurons [19,20]. When the moment method 
is applied to a single neuron model with K variables, X-dimensional stochastic DEs are 
replaced by {1/2)K(K + 3)-dimensional deterministic DEs. When the moment method is 
applied to A^-unit neuron ensembles under consideration, A'A^-dimensional stochastic DEs 
are replaced by Agy-dimcnsional deterministic DEs where N^q = {1/2)KN{KN + 3) [17]. 
For example, in the case of AT = 2 (FN model), the number of equations is Ngq = 230, 20 300 
and 2 003 000 for N = 10, 100 and 1000, respectively. In the case of AT = 4 (HH model), we 
get Neg = 860, 80 600 and 8 006 000 for N = 10, 100 and 1000, respectively. These figures 
are too large for us to make simulations for realistic neuron clusters. In their subsequent 
paper [19], RT transplanted the result of the moment method for HH neuron ensembles to 
FPE- type equation which has not been solved yet. 

In a previous study (Ref. [21] is hereafter referred to as I), the present author proposed a 
semi-analytical dynamical mean-field approximation (DMA), in which equations of motions 
for means, variances and covariances of local and global variables were derived for A^-unit 
FN neuron ensemble. The original 2A^-dimensional stochastic DEs are replaced by eight- 
dimensional deterministic DEs: N^q = 8 is much smaller than corresponding figures in the 
moment method mentioned above. DMA calculations in I on the spiking-time precision and 
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the synchronization in FN neuron ensembles are in good agreement with direct simulations. 
The feasibihty of DMA has been demonstrated in I. 

The purpose of the present paper is two folds. The first purpose is to extend DMA of I to 
general neuron ensembles subject to white noises described by XA^-dimensional stochastic 
DEs, which will be replaced by K{K +2)-dimensional deterministic DEs. The second purpose 
of the present paper is to apply the generalized DMA to an ensemble of HH neurons, which 
is more realistic than FN neuron model previously studied in I. Since Hodgkin and Huxley 
proposed the HH model in 1952 [1], much studies have been intensively made on properties 
of the HH model. Responses of a single, pairs and ensembles HH neurons mostly to direct 
and sinusoidal currents have been investigated. In recent years, responses of HH neurons to 
spike-train inputs have been studied [22]- [25]. The stochastic resonance (SR) of HH neurons 
for sinusoidal and spike inputs with various kinds of added noises has been investigated [26]- 
[33] . These studies have shown that noise can play a constructive role in signal transmission 
against our conventional wisdom. In most studies on SR, however, noises added to ensemble 
neurons are considered to be independent of each other. Quite recently effects of spatially 
correlated noises on SR have been investigated [31], which shows that although common 
noises work to enhance the synchronization in neuron ensembles, they are not effective for 
SR, in contrast to independent noises. We will adopt in this study, spatially correlated 
white noises in order to clarify respective effects of common and independent noises on the 
response of ensemble neurons. 

The paper is organized as follows: In Sec. II, we extend a DMA theory to general 
neuron ensembles described by KN stochastic DEs. Our DMA theory is applied to HH 
neuron ensembles in Sec. HI. Some numerical results on HH neuron ensembles are presented 
in Sec. IV. Conclusions and discussions are given in Sec. V. 

II. DMA FOR A GENERAL NEURON ENSEMBLE 

A. Equation of motions 

We assume an ensemble of A^-unit neurons (A^ > 2), each of which is described by K- 
dimensional non- linear differential equations (DEs). Dynamics of a given neuron ensemble 
is expressed by 

^ = F^'\M) + (^) E G{vM + K^'\t)+Ut), (1) 

^-F(^\M), (p = 2toX) (2) 

where Vi = Upi with p = 1 denotes the membrane potential of a neuron i (= 1 to A^), 
Upi with p = 2 to K stands for auxiliary variables and F^^^ is functions of {{uqi}). The 
synaptic-coupling strength w is assumed to be constant, G{v) = 1/[1 + exp[— (f — 9)/e]] is 
the sigmoid function with the threshold 9 and the width e [34,35], and K^^^ stands for an 
applied external input whose explicit form will be given later [Eq. (56)]. The last term of 
Eq. (1) expresses the spatially correlated white noises, ^i{t), given by 
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< m > = 0, (3) 

^{I3l + (3j5,^)5{t-t'), (4) 

where j3c = Pi and (3i = \J Po ~ Pi denote the magnitudes of common and independent 
noises, respectively, and the bracket < • > expresses the stochastic average [36]; the case of 
Pi — {Pi — Po) stands for independent (common) noises only. 

In order to derive DEs in DMA theory, we first define the global variables for the ensemble 
by [21] 

U,{t) = (l/N) Y: (5) 

i 

and their averages by 

Mt)=iiu,it)=<Up{t)> . (6) 
Deviations from these averages of local variables are given by 

6upi{t) = Upiit) - i^up{t), (7) 
and those of global variables given by 

6Up{t) = Upit)-fiuM- (8) 

Next wc define the variances and co- variances between local variables given by (argument 
t is neglected hereafter) 

i 

and those between global variables given by 

Pp,q = Pup,Ug =< 5Up SUq >, (10) 

It is noted that 7«p,«q expresses fluctuations in local variables while Pup,uq those in global 
variables. 

We assume that the noise intensity is weak and that the distribution function p{z) for 
i^A^-dimensional random variables of z = ({wpj}) is given by the Gaussian distribution 
concentrated near the mean point of /i=({/i„p}) [36]. Numerical simulations have shown 
that for weak noises, the distribution of v{t) of the membrane potential of a single HH 
neuron nearly obeys the Gaussian distribution, although for strong noises, the distribution 
of v{t) deviates from the Gaussian, taking a bimodal form [22] [37]. Similar behavior of 
the membrane-potential distribution has been reported also in a FN neuron model [18] [38]. 
By using Eq. (7), we express Eqs. (1) and (2) in a Taylor expansion of 6upi up to the 
fourth-order terms. The average yields DEs for the means of djiu^/dt [Eq. (16)]. DEs of 
variances and covariances may be obtained by using the equations of motions of 5upi. For 
example, DE for djup,uq/dt is given by 
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with 

q q r ^ q r s 

+ V Hi + E + \Gu,uAKk - 7i,i) + \Gu,u,u, Kk\ + ^^'^}, (12) 

where g, r and s run from 1 to K, F^"'') = F^^f) = dF^P^/dug, F^^^^l = d"^ F^p^ du^dur 
and F^^SrUs ~ F'^py dugdurdus evaluated at the means of {{f^up}), and similar derivatives 
for G. In the process of calculations of means, variances and covariances, we have taken into 
account the fourth-order moment contributions with the use of the Gaussian approximation, 
as given by 

< 6Upi6Uqi6UriSUsi > ~ < 6Upi6Ugi >< dUriSUsi > + < SUpi6Uri >< 6Uqi6Usi > 

+ < SUpiSUsi >< SUqiSUri >, (13) 

i 

E E ^ SUpiSUqjSUrjSUgj > — Pup,Uqlur,Us Pup,Urluq,Us + Pup ,Us 7wq ,«r- ) (l^) 

« i 

The importance of including the fourth-order term has been pointed out by Tanabe and 
Pakdaman [38] in the improved moment method for a noisy FN neuron. 

After some manipulations, we get DEs for means, variances and covariances given by 
(details being given in Appendix A of I): 



^^-F(-'^+li:i:Fi::huq,ur+Sp, [wUo+k(% (i6) 



q r 



dt 



r 



r s w 



dpup, 



/ J I 11, • r u,q,u,r ' - rup,u,ri ' L^/'U ' y jy 

+W Ul[5pl Puq,Ul + Sql Pup,Ul\ 

E E '^2lFi^ulutiP'^q'UrlfUs,Ut + Puq,UslUr,Ut + Puq,Utlur ,Us) 

^ r s t 

'^F^^!l}^y^^{pUp,Urlus,Ut + Pup,Uslur,Ut + Pup ,Ut Tu,- ,Us ) ] ) (l^) 



with 
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^0 = ^ E < G{vj) >= G + ^G,,7.,., (20) 

j 

Ul — Gy + -^VVV1V,V1 (21) 

where expresses output spikes of the ensemble, = uy, and arguments of r, s and w in 
the sums run from 1 to K. The original K A'"-dimensional stochastic DEs are transformed 
to A^eg-dimensional deterministic DEs where N^g ^ K + K{K + 1) = K{K + 2). 



B. Property of DMA 

In previous Sec. IIA, DMA has been derived with the use of equations of motions 
for moments. It is, however, possible to alternatively derive DMA from the conventional 
moment method with a reduction in numbers of variables, as was shown in I for FN neuron 
ensembles. In Appendix A, we present a derivation of DMA from the moment method for 
a general neuron ensemble under consideration. 

We should note that the noise contribution is /9g in Eq. (17) while that is [(l/A^)/9g + 
(1 — l/N)Pf] in Eq.(18). When model parameters of Po, Pi, w and are varied, the ratio 
of Pv,v/lv,v changes. In particular, in the case oi w = 0, we get 

(22) 

(23) 
(24) 

Equation (23) agrees with the central- limit theorem for independent noises while Eq. (24) 
expresses the result for common noises. On the other hand, in the opposite limit of w — > oo, 
we get Pv,v/lv,v 1- The change in the ratio of Pv,v/lv,v reflects on the the flring time 
distributions and the degree of synchronization in neuron ensembles, as will be discussed in 
the foUowings. 

Firing Time Distributions 

The nth firing time of a given neuron i in the ensemble is defined as the time when the 
membrane potential Vi{t) crosses the threshold 9 from below: 

toin = {t\vi{t) = e;Vi>Q}. (25) 

The distribution of firing times of tom a given neuron i is given by [17] [21] 

m-<f>C-^) i(^) ©(/i.), (26) 




with the normal distribution function given by 
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(f){x) 



1 / 



and 



27r 



(27) 



(28) 



where cr^ = ^^/'Jv,v and ^t^, = dfi^/dt evaluated at t — where i^vitl) — 0. In the hmit of 
vanishing Eq. (26) reduces to 



z,{t)^5{t-t:). 



(29) 



Similarly we may define the mth firing time relevant to the global variable V{i) 
{l/N) EiViit) as [21] 



tgm = {t\V{t) = 9;V{t)>0}. 
The distribution of firing times of tgm is given by 

f t — t*o\ d f f^i 



with 



5tog J dt 



(30) 



(31) 



(32) 



where ag — y/Pv,v In particular, in the case of no couplings, we get 



— - l + fl-l) 



(for w — 0) 



(33) 



Synchronous Response 

The synchronization ratio is defined by [21] 



(1 - l/N) 7^,^ 



with 



Cv,v = i^Pv,v - lv,v) = < >' (35) 

expressing the averaged covariance for the variable of {{5vi}). S{t) changes as the model 
parameters of /3o, f3i, w and are varied. It is easy to see from Eqs. (23) and (24) that 
S = (the asynchronous state) for w = and Pi -C (3o, while S = 1 (the completely 
synchronous state) for w ^ /?q or /3i = (3q. In particular, for w = 0, we get 

S{t) = (^^^ , for ^ = (36) 

which implies that the synchronization is induced by common noises. 
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III. DMA FOR HH NEURON ENSEMBLES 



A. Equation of Motions 

For the HH neuron model (i^ = 4), F^^) in Eq. (1) is given by [1] [23] 

^(1) ^ F^^Xvi, rrii, hi, rii) = -^[gNa.m^hi{vi - ^Na) + QKntivi - vk) + Qhivi - vl)], (37) 
jrip) ^ Upi) = -K^ivi) + b^^ivi)] u^^ + a^^{vi), (p = 2 to 4) (38) 

In Eqs. (37) and (38) Uu = Vi expresses the membrane potential of a neuron i, and U2i = rrii, 
Usi — hi and u^i = rii denote gate variables of Na and K channels for which aup{v) and bu^{v) 
(p = 2 to 4) are given by 

_ 0.1 (^; + 40) 
""^v^^ ~ [1 _ g-(^;+4o)/io] ' v39) 

6„(^) = 4e-(^+^^)/^^ (40) 
a^(^;) = 0.07e-(''+*^5)/'o, (41) 



^'^(^) = [1 + e-(^+35)/io] • (42) 
0.01 (v + 55) 

«-(^) = [l_e-(^+55)/io] ' (43) 

6„(^;) = 0.125 e-(''+^5)/«°. (44) 

In Eq. (37), the reversal potentials of Na, K channels and leakage are v^^ = 50 mV, 
Vk = —77 mV and = —54.5 mV: the maximum values of corresponding conductances 
are g^^ = 120 mS/cm^, qk = 36 mS/cm^ and gi^ = 0.3 mS/cm^: the capacitance of the 
membrane is C = 1 /xF/cm^. From functional forms for F^'"^ and F^'^p^ given by Eqs. (37)- 
(44), we get F,,, = 0,Fi«^) = F^^^) S.^F^^^^ = F^^' S,, and F^^ = 0. Numbers of 

non- vanishing third-order derivatives are six for F^'') [F^^^, F!;^2h, F^±, F^l^, F^± and 
Fit,] and two for each F(«-) (p = 2 to 4) [Fi«„-) and F,t^)]. 

After some manipulations with Eqs. (16)-(18), we get DEs for means, variances and 
covariances given by (p, g = 2 to 4) 

^ - F^'' + ^ E E Fili^ lu,,u. + E ^i£7">«. + ^ ^0 + K^'\ (45) 

^ p=2 q=2 p=2 



8 



(it „ o 

9=2 



Up ^ Uq )lUp,Uq^^v IV,Uq~^v 1 V ,Up ^ ^^Up,Uq1 V*^ ) 

2[i^i"V.,. + E Fi^^ P.,up] + 01 + (1 - ^)/3i'] + 2 ^ C/i p.,. + (50) 



.^^(^,) iV-pM^ + F^^f^O +wU^O +Y (51) 

1, K-'^V ' ^Up ) Hv,Up ^ J^Uq HUq,Up ^ V HV,V T W U \ Pv,Up ^ ^V,Up, \^^) 

q=2 

^P'^p'^i - (riup) , pK)^^ _L p(«p)^ _L +Y (52) 

~ Up ^ Uq )h'Up,Uq \ y HV,Uq \ v HV,Up ^ -< Up,Uq) ) 



with 

Cup.Uq — ^ — 1 ) P'<J'P,1J-q ~ lup,Uq)-l (53) 

^0 = ;^ E < ^(^^O >= G + ^G,,7,,„ (54) 
i 

^^1 = + -G„uu7^,^, (55) 

where F'^'"\ F^'"^ — dF^'"'> /dv et al. evaluated at means of {iJ,y, /j^m, Ph, Pn)- In Eqs. (45)-(52), 
Xy^y and Yy^y et al. denote the contributions from the fourth-order terms, whose exphcit 
expressions are given by Eqs. (B1)-(B6) in Appendix B bccuase they are rather lengthy. 
Although calculations of the fourth-order terms are rather tedious, they play important roles 
in stabihzing DEs. This is numerically demonstrated in Appendix B for the case of = 1. 

The original 4A"-dimensional stochastic DEs given by Eqs. (37) and (38) are transformed 
to 24-dimensional deterministic DEs given by Eqs. (45)-(52) with Eqs. (B1)-(B6): four 
means, (p^,, /i™, /i/t, Ai„), ten moments for local variables (7^,,^,, 7™,™, 'yh,h, ln,n, lv,ra, lv,h, 
lv,n, lm,h, 7/i,n, 7m,n)) and tcu momcuts for global variables (p^,^, Pm,m, Ph,h, Pn,n, Pv,m, Pv,h, 
Pv,ni Pm,,hi Ph,ni Pm,n) ■ 

In this subsection, DMA for HH model has been obtained by the method of equations of 

motions of means, variances and covariances of local and global variables. Wc may, however, 
derive it from the moment method, as mentioned before. In Appendix C, DEs in the moment 
method are presented for HH model. 

We expect that our DMA equations given by Eqs.(45)-(52) and (B1)-(B6) may show 
much variety depending on model parameters such as the strength of white noise (/3o, Pi), 
couplings (w) and the ensemble size (N). In the next Sec. IV, wc will present some numerical 
DMA calculations, which are compared with simulation results. DMA equations have been 
solved by the fourth-order Runge-Kutta method with a time step of 0.01 ms for the initial 
conditions of p„ = -65.0, p^ = 0.0528, Hh = 0.597, pn = 0.317, and 'yup,uq = Pup,uq = 
{up, Uq=v, m, h and n). Direct simulations have been performed by solving 4A"-dimensional 
DEs given by Eqs. (37) and (38) by using also the fourth-order Runge-Kutta method with a 
time step of 0.01. Simulation results are the average of 100 trials otherwise noticed. 
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IV. CALCULATED RESULTS OF HH NEURON ENSEMBLES 



A. Firing time distribution 



In the present study, we pay our attention to the response of the HH neuron ensembles 
to a single spike input applied to all neurons in the ensemble, given by 



where 0(,t) = 1 for a; > and otherwise, stands for the magnitude of an input spike, 
C the membrane capacitance [Eq. (37)], ti the input time of a spike, and (=1 ms) the 
time constant of synapses. We get the critical magnitude of Jjc = 3.62 /xA/cm^, below which 
firings of neuron cannot take place without noises (/3o — Pi — 0). We have adopted the value 
of li —5 ^K/cFo? for a study of the response to a supra-threshold input. We express the 
couphng constant why w = J /C with J in units of jiA/cm^. The time, voltage, current and 
noise intensity are hereafter expressed in units of ms, mV, /lA/cm^ and V/s, respectively, 
though they are sometimes omitted for a simplicity of our explanation. We have adopted 
parameters of ^ = mV and e = 10 mV in the sigmoid function G{v) such that output Uq is 
similar to the result given by the alpha function [see Fig. 1(a)]. Adopted parameter values 
of /9o, J and will be explained shortly. 

Figures 1(a), 1(b) and 1(c) show the time courses of /i^, (= ^J^|v,v) and ag (= y/Pv,v), 
respectively, when a single spike is applied at t = 100 ms. Solid and dashed curves express the 
results of DMA and direct simulations, respectively, which are calculated with parameters 
of /3o = 0.1, /3i = 0, J — and N — 100. States of neurons in an ensemble when an input 
spike is injected at t = 100 ms, are randomized because noises have been already added since 
t = 0. We note that fi^ obtained by DMA is in very good agreement with that obtained by 
simulations as shown in Fig. 1(a), where an external input of K^^^t) and an output of Uo{t) 
are also plotted. Figures 1(b) and 1(c) show that and ag calculated by DMA are again 
in good agreement with those of simulations. We note that the relation given by Eq.(22): 

= CTi/N valid ior w = J/C = 0, is supported by our calculations. 

Figure 2(a) shows Zi, the firing probability of local variable, which is calculated for 
Po = 0.1, Pi = 0, J = and = 100. Firings occur at t ~ 103.6 ms with a delay of about 
3.6 ms. Fluctuations of firing times of local variable, 5toi, are 0.066 ms in DMA while it is 
0.069 ms in simulations which is the root-mean-square (RMS) value of firing times defined 
by Eq. (25). In contrast. Fig. 2(b) shows Zg, the firing probability of global variables. 
Fluctuations of firing times of global variable 6tog are 0.0066 ms in DMA and it is 0.0083 
ms in simulations, respectively. We note that 5tog is much smaller than 6toi [Eq. (33)]. 

Noise-strength dependence 

When the noise strength is increased, the distribution of membrane potentials is widen 
and fiuctuations of firing times are increased, as was discussed in Sec. IIB. Filled squares in 




(56) 



with the alpha function: 




(57) 
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Fig. 3(a) show the (3o dependence of Stoi obtained by DMA theory with /?i = 0, J = and 
N = 100, while open squares express the RMS value of firing times obtained by simulations. 
The agreement between the two methods is in fairly good for Pq < 0.1 but becomes worse 
for P > 0.1. In contrast, filled circles in Fig. 3(a) show the /3o dependence of Stgg relevant to 
the global variable obtained by DMA theory and open circles stand for RMS values of firing 
times in simulations. We note that Stag is much smaller than Stoe because Stag — Stge/VN 
[Eq. (33)]. 

As Pi is increased for a fixed Pq, the contribution from common noises increases while 

that from independent noises decreases {Pc = Pi, Pi = \j Pi — Pi)- The Pi dependence 
of firing-time fiuctuations is shown in Fig. 3(b). Filled squares and circles denote the 
results of tog. and tog, respectively, obtained by DMA, and open squares and circles those 
by simulations. Figure 3(b) shows that 5tog is almost linearly increased as Pi is increased, 
while Stoi, remains constant. In the limit of Pi = Po = 0.1, for which only common noises 
are applied {Pc = 0.1 and /?/ = 0), we get Stag = Stoe, which shows that common noises do 
not work to reduce global fiuctuations. 

Ensemble-size dependence 

Filled squares in Fig. 4(a) show the N dependence of Stoe relevant to local fluctuations 
for Po = 0.1, Pi = and J = 0, obtained by DMA theory, while open squares express that 
obtained by simulations. We note that Stoe is independent of N because of no couplings 
{J — 0). In contrast, Stog relevant to global fiuctuations inversely decreases when the size 
is increased, as shown by filled and open circles which are obtained by DMA theory and 
simulations, respectively. The relation: Stog oc holds as given by Eq. (33) for 

Pi = 0. Figure 4(b) shows a similar plot for a finite value of Pi = 0.05 with Pq = 0.1, J = 
and A" = 100. In the limit of A^ — > oo, the ratio of Stog /Stoe approaches a finite value of 
Pi/Po^O.5 [Eq. (33)]. 

Coupling-strength dependence 

So far we have neglected the coupling of J, which is now introduced. Filled squares in 
Fig. 5(a) show the J dependence of Stot calculated by DMA theory for Pq = 0.1, Pi ~ and 
A^ = 100, while open squares that obtained by simulations. Filled and open circles express 
Stog in the DMA theory and simulations, respectively. We note that Stoe is much reduced as 
J is increased although there is little change in Stog. Figure 5(b) shows a similar plot of the 
J dependence of firing time accuracy for finite Pi = 0.05 with Po = 0.1 and A^ = 100. Again 
a reduction in Stoe ^ increasing J is more significant than that of Sog. 

B. Synchronization ratio 

One of important effects of the couplings is to yield synchronous firings in ensemble 
neurons. Figures 6(a) and 6(b) show the time course of the synchronization ratio S{t) for 
J = 100 and 200 /xA/cm^, respectively, with Po = 0.1, Pi = and A^ = 100: solid and dashed 
curves denote the results of DMA and simulations, respectively. Fairly large fiuctuations in 
simulation results are due to a lack of trial number of one hundred, which is a limit of our 
computer facility. A comparison between Figs. 6(a) and 6(b) shows that S{t) is increased 
as J is increased: the maximum value of S{t) in Fig. 6(b) is Smax — 0.019 which is larger 
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than Srnax = 0.007 in Fig. 6(b). Figure 6(c) shows the time course of S{t) for a finite 
Pi = 0.05 with /3o = 0.1, J = 100 and N = 100. A significant increase in S is reahzed 
at 100 ~ t ~ 120 ms which is induced by an apphed spike [note the difference in vertical 
scales of Figs. 6(a), 6(b) and 6(c)]. We note a fairly large value of 5" = 0.25 even without 
an applied input spike at i ~ 100 or i ~ 120. This expresses the synchronization among the 
membrane potentials of ensemble neurons induced by added noises although they do not 
induce firings. In order to distinguish the synchronization with firings from that without 
firings, we define the firing-induced synchronization ratio, S'{t), given by 

S'it) = Sit) - S,, (58) 

where Sb = {Pi/PoY denotes the background synchronization induced by noises only [Eq. 
(36)]. We get Smax — 0.369, S'^^^^ — 0.119 and Si, = 0.25 in Fig. 6(c). From a comparison 
of Fig. 6(c) with Fig. 6(a), we note that S'{t) is also much increased by common noises. 

An increase in S{t) by an increase of J is clearly shown in Fig. 7(a), where the maximum 
of S{t) (Smax) is plotted as a function of J. A disagreement between results of DMA and 
simulations for J < 50 is due to fluctuations in simulations because of insufficient trial 
number as mentioned above. The dependence of S^ax on the size N is shown in Fig. 7(b) 
where /3o = 0.1, /3i = and J = 100. S^ax is decreased as increasing N. Figure 7(c) expresses 
the po dependence of S^ax for pi = 0.05, J = 100 and N = 100. At pQ = pi = 0.05, we get 
Smax = 1, which is decreased as increasing Pq. Filled squares in Fig. 7( c) denote S'^^^^^, which 
shows the maximum around Pq ~ 0.08. In contrast, Fig. 7(d) show the Pi dependence of 
Smax for /^o = 0.1, J = 100 and N = 100. Smax is increased as increasing Pi, and approaches 
unity as Pi — > po {— 0.1). We note that S'^^^^^ has the maximum at Pi ~ 0.07. 



V. CONCLUSION AND DISCUSSION 



In previous Sec. IV, we have reported DMA calculations for a single spike input to HH 
ncmon ensembles. DMA calculations and simulations have shown that 

(a) 6to£ increases as increasing Po, or decreasing J, independently of Pi and N, 

(b) Stog increases as increasing Pq or Pi, or decreasing N, independently J, and 

(c) Smax increases as increasing Pi or J, or decreasing Pq ot N. 

In order to understand these behaviors, we have tried to obtain phenomenological, ana- 
lytical expressions for 6to£, St^y and Smax as functions of Po, Pi, J and N. For small J, we 
express 7t,,„ and p^,^ in power series of J at t = t* where neurons fire, given by (see Appendix 
E of I, [34]) 



lv,v (x Po 1- (aiJ + a2 + ••) 



Pv,v OC Po 



N 



+ 1 




(59) 
(60) 



It is noted that in the hmit of J = (w = 0), Eqs. (59) and (60) reduce to Eq. (22). 
Substituting Eqs. (59) and (60) to Eqs. (28) and (32), we get 
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Stoe oc Po 

Stag OC Po 



l-^(aiJ + a2J^ + --) 




1/2 



(61) 
(62) 



Equations (61) and (62) may explain the behavior of Stoe and Stog in items (a) and (b) 
mentioned above. 

Next we will obtain the analytical expression for Smax- For small J, we get 7^^^ and p^^y 
m power series of J at t = 4""^ where S{t) takes the maximum value, given by (see Appendix 
E of I, [34]) 



7^,„ oc /^o^ [1 - (biJ + b2J^ + ■■) 
Pv,v oc /?o ^ + ( 1 




N 



Substituting Eqs. (63) and (64) to Eq. (34), we get 



+ 



(63) 
(64) 





hJ+{b2 + bl)J^ + 



(65) 



Equation (65) is consistent with the item (c): Smax increases as increasing f3i and J, or 
decreasing Pq and N. In the case of j3i = 0, Eq. (65) shows that Smax is independent of j3o, 
which is supported in DMA calculation and simulations (not shown). Expressions given by 
Eqs. (61), (62) and (65) are useful in a phenomenological sense. In principle, expressions as 
given by Eqs. (59) and (60) may be derived from DMA equations given by Eqs. (45)-(52) 
although we have not unfortunately succeeded in getting them because of their complexity. 

Numerical calculations in Sec. IV have been made for the response to a single spike 
input. DMA is, however, applicable to arbitrary inputs. This will be demonstrated by 
adding spike trains to HH neuron ensembles, given by 



(66) 



where tj„ expresses the nth input time. Figures 8(a) and 8(b) show the time courses of fi^ and 
{=^yiv,v), respectively, for Poisson spike trains with the average interspike interval (ISI) 
of 25 ms; solid and dashed curves express results DMA and simulations, respectively, for 
/?o = 0.1, /3i = 0, J = and N = 100. The time course of /ly of DMA is in good agreement 
with that of simulations. A comparison between the input K^^'' and output Uq shows that 
when ISI of input is shorter than about 10 ms, HH neurons cannot respond because of the 
refractory period. Figures 8(b) shows that of DMA is also in good agreement with that 
of simulations. 

To summarize, the DMA theory previously proposed for FN neuron ensemble in I, has 
been generalized to an ensemble described by XA^-dimensional stochastic DEs, which has 
been be replaced by K(K + 2)-dimcnsional deterministic DEs expressed by means and 
second-order moments: contributions from the fourth-order moments are taken in account 
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by the Gaussian decoupling approximation. DMA has been apphed to HH neuron ensem- 
bles, for which we get 24-dimensional deterministic DEs. We have studied effects of noise, 
the coupling strength and the ensemble size on the firing time precision and the firing syn- 
chronization for single-spike inputs, obtaining the following results: 

(i) the firing-time accuracy of the order of one-tenth ms is possible in a large-scale HH 
neuron ensemble, even without couplings, 

(ii) the spike transmission is improved with the synchronous response by increasing the 
coupling strength, and 

(iii) the synchronization is increased by common noises but decreased by independent 
noises. 

The item (i) and (ii) are consistent with the SR results in HH neuron ensembles [28]- [33]. 
Although they are quite similar to the case of FN model discussed in I, their quantitative 
discussions are possible with the use of the realistic HH model. The item (iii) agrees with 
the result of Ref. [31] for SR in HH neuron ensembles subject of common and independent 
noises. 

Our calculations have demonstrated the feasibility of DMA, whose advantages may be 
summarized as follows: 

(1) because of the semi-analytical nature of DMA, some results may be derived without 
numerical calculations, 

(2) DMA is free from the weak-coupling constraint although it assumes weak noises, 

(3) a tractable small number of DEs makes calculations feasible for large-scale neuron 
ensembles with a fairly short computational time, 

(4) DMA may be applicable to ensembles with fluctuations not only due to noises but 
also due to some inhomogeneities in model parameters, and 

(5) DMA can be apphed to more general stochastic systems besides neuron models. 

As for the item (3), we may point out that for example, the CPU time of DMA calcula- 
tions for a 200 ms time course of a = 100 HH neuron ensemble with the use of 1.8 GHz 
PC is 2 s, which is about 2500 times faster than the CPU time of 85 min (~ 5000 s) for 
direct simulations with 100 trials. It is necessary to stress the importance of the fourth- 
order contributions in stabilizing solutions of DMA, which is numerically demonstrated in 
Appendix B. Although expressions for the fourth-order contributions are lengthy, we are 
much beneflted from them once they are derived and planted into computer programs [39]. 
This paper is the flrst report on our DMA calculations of HH neuron ensembles. We are now 
under consideration to incorporate the time delay in the coupling terms in Eq. (1), with 
which the HH neuron ensemble may show the intrigue behavior like chaos. Such calculations 
will be reported in a separate paper. 
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APPENDIX A: DERIVATION OF DMA FROM THE MOMENT METHOD FOR 

GENERAL NEURON ENSEMBLES 



In the moment method, we define the means, variances and covariances given by [17] 

Cit^. = < A«p.A«,,- >, (A2) 

where Aitpj = Upi—m\^. Assuming the weak couphngs and adopting the Gaussian decouphng 
approximations for the fourth-order moments, we get DEs for general neuron ensembles 
described by Eqs. (1) and (2): 



dt 



i:iG + lG^,,u,,CtiJ + K^% (A3) 



1^ Y^ Y^ Y^rci(tjp) _i_ f<i,i f<i,j , f<i,i fij \ 

Q ' ' / / '\ UriUsiU^i\^Ur,Us Uiu,Uq ~ ^Ur,Uw^Us,Uq ~ ^Us,Uw Ur,Uq/ 



r s w 



where F^"^) = = dF^P^dUri and F^pJ.^^^^^. = 9(3)F(^')/aMH<9M«<9M«,i evaluted for 

the means of ({m^ }), and the last term in Eq. (A4) denotes the fourth-order contribution. 
The number of DEs is N^^ = KN + {1/2)KN{KN + 1) = {1/2)KN{KN + 3). 

In order to derive DMA from the moment method, we define the quantities given by 



N 

%,a = ^E^k:\+<a, (A6) 

Pk,x = 4^ > {K,X^Up, Ug) ( A7) 
j 



where 



<^ = J^Y.^<^<^ (A8) 

i 

5m\ ^m\-n^. (A9) 

We may show that Eqs. (A3) and (A4) with Eqs. (A3)-(A7) yield Eqs. (16)-(18): Jl^ = 
7k,a — 7k,a ^iid Pk,a — Pk,a- Then the moment method yields the same results as DMA as 
far as the averaged quantities are concerned (see also Appendix B of I) . 
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APPENDIX B: THE FOURTH-ORDER CONTRIBUTIONS IN DMA FOR HH 

NEURON ENSEMBLES 



The fourth-order contributions given by X^^x and Y^^x {i^,^ = v,Up) in Eqs. (45)-(52) 
are expressed by 

-^v,v Pvrnm^1v,v1m,m ~\~ '^Ivymlvj-m) ~l~ Fvmh^'~1v,v'~1m,h ~\~ '^1 v ,m'~1 v ,h) 

,v'yn,n ~l~ '^'yv,n'yv,n) -frreryirTiTf ,mTm,m ~l~ Pnnn1v,n^n,n 
+ F!^lhhv,hlm,m + '2lv,mlm,h) , (Bl) 

-^v,Up 2^-^^mmi'y'v,Up'ym,m ~l~ '^lfup,mlfv,m) ~l~ ■^vmh^'^v,Up1m,h ~l~ lfup,mlfv,h ~l~ lfup,hlfv,m) 

'^■^vnni'JvjUpJn.n + '^^Up,n1v,n) ~l~ F^'^^^^^^ ^i'^jf^.^y^ + F^^^'~)up,ii)n,n 

+Ftal ilv,uplv,v + 2p.,.7.,«p)], (B2) 

^,11 — -^imm (P'fifTmjm + '^Pv,m')v,rn) + F^jj^f^{pi,^v^m,h + Pv,mlfv,h + Pv,h1v,m) 

+^mL(p^,/»7m,m + 2p^,r„7m,/i), (B4) 

^,Up = 2 '^Pup,m'Jv,m) + F^j^j^[py^up1im,h + Pup,m1v,h + Pup,h1iv,m) 

~'~-^inn(P'!',«p7n,n + '^Pup,n1v,n) + -f^m,rnm,PMp,"i7m,m + -^nr!.nP«p,n7n,n 
'^F^j^li{^pUp,h1m,m + ^Pup,m'ym,h) + F^^^^ p^^^^y^^ 

)], (B5) 

Yup,Uq 2 ['^'^^'^^''^'"''"9^'^''^ ^ Fyy[J^{Pup,Uqlfv,V ~l~ Pv ,Uqlv ,Up^ ~^ F^^!^^ Py^Up'~1v,V 

'^Fy^!^^{pUp,Uq1v,V + '^Pv,Uplv,Uq )] (B6) 

where F^^^/j = d^F^'"^ / dvdmdh et al. Although calculations and computer programming of 
fourth-order contributions given by Eqs. (B1)-(B6) are rather tedious, they play important 
roles in stabilizing the solution of DEs [39]. 

Here we demonstrate the importance of the fourth-order contributions in the case of a 
single HH neuron (iV = 1) for which w = 0, 7«;,a = Pk,x and = ^k.a in Eqs. (45)-(52) 
and (B1)-(B6). Figure 9(a) shows the time course of p^ for I3q — 0.1 and /3i = when the 
constant input of = 10 /lA/cm^ is applied at t = ms. The solid and dashed curves 
express the results of DMA and the simulation (100 trials), respectively. The dotted curve 
denote the result of DMA2 (the second-order DMA) in which the fourth-order contributions 
are neglected {X^^x = Yf^^x = 0). For t < 60 ms, all results seem to be in good agreement. 
At i ~ 60, however, the solution of DMA2 becomes unstable and significantly deviates from 
those of DMA and the simulation. From the time course of ae — shown in Fig. 9(b), 

we note that such deviation of DMA2 ah^eady starts from t ~ 30 ms. The solution of DMA2 
is stable at /? < 0.037 for the constant current of /j = 10 pK/cw?. 



16 



Figures 10(a) and 10(b) show the time courses of ^.y and for (3q = 0.2 and (3i = when 
we apply the periodic spike train input given by Eq. (66) with li = 5 jiA/cw? and a constant 
ISI of 25 ms. Figure 10(b) clearly shows that the result of DMA2 deviates from those of 
DMA and simulations from the first spike input and that the result of DMA2 diverges at 
the second spike input. The solution of DMA2 is stable only at /3 < 0.178 for this periodic 
spike. 



APPENDIX C: THE MOMENT METHOD FOR HH NEURON ENSEMBLES 

We will derive DEs in the moment method for HH neuron ensembles, defining the means, 
variances and covariances given by [17] 

w}^^<Vi >, (CI) 

^<y'pi>, (C2) 

Qi = < AviAvj >, (C3) 

= < Av.Aup, >, (C4) 

Cii^, = < ^^pAu,j >, (C5) 

where Avi = Vi — ml and Aupi = Upi — m\^^. Adopting the weak-coupling approximation 
and the Gaussian approximation for the fourth-order terms, we get DEs given by 



-t 4 4 'i 
_|_ _ C^<^ _|_ pMr<i,i 

Q / ^ / ^ Up,Uq Up,Uq ~ / y V,U V,Up 

^ p=2 q=2 p=2 

+ (]^) E iG^""'^ + iGt'^C'^:^) + K^'\t), (C6) 



dml 



+ lir(%>Oc:;; + (C7) 



dt 2 



dt 



dt 



Up ,Uq 



dt 

with 



p=2 

+ [ E G^:"^ Ciit + E G^:''^ + (C8) 

q=2 

+ (]^) E G^:^^ + , (C9) 

+ ) + + F^'^^'^Cii + Zii\ , (CIO) 

\ Ui-p U/q ^ Up U/q (JLp J U/q / U Ctq (/ f , Up Up ^ Ug ' \ / 



— 1 [ fr(.v) (fiJ fjJ I OC*'-' C-''-' -I- C*'-^ C*'* -I- 2r'-''* C*'-' 

Q L vmm\ v,v m,m ~ ^v,m v,m ~ v,v m,m ~ ^v,m v,r, 
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' vmh\ v,v m,h ' v,m v,h ' v,h v,m ' v,v m,h ' v,m v,h ' v,h v,mJ 
~ vnn\ v,v n,n ~ v,n v,n ~ v^v n^n ~ v,n v,n) 

-lpC^) f/^^J /^JJ -L r^*'* -L pC^) (n'i'd r<3,j \ r<j,i r<i,i \ 

I f(^) (f^iJ f^jj lOC^'^C-^'^ \ fhi r<i,i _|_ or"-?'* r**'* W (TlT) 
' rnrnh\ v,h rn,rn ' v,m m,h ~ v,h rn,rn ' v,m rn,h)\^ \^ J 

r7i,j _\\-p{v) ni,i I 2C^'j _|_ (C'^d n^'^ + (7*'-^ (7*'* +(7*'-' (7*'* 1 

v,Up c^V vmm\ v,Up m,m ~ m,Up v,mJ ~ vmh\ v,Up m,h ~ m,Up v,h ~ h,Up v,m) 

_l_ Tpi-") (nhj + 2<7*'-' \ + p(^) <7*'-' <7*'* + /^(^) <7*'^ <7*'' 

' t)nnV^ii,Up n,n ' n,Up v,nJ ' mmm m,Up m,m ' nnn^n,Up^n,n 
' mmh\ h,Up^ m,m ' ^rn,Up^m,h) ~ tJiiii 

+^2i(C::i,C^;i + 2C:;iC^;iJ], (C12) 

Up,Uq 2^ I VVV ^V,Uq V,V ' VVUp\ Up,Uq V,V ' V,Uq V,Up) ~ VVV Up,V V,V 

+F^iiCi;l^^Ci;^ + 2C^,„C^;iJ], (C13) 

where F^'"^\ F^^^^ = dF/dvi and -Fj'^^ = d'^F/dvidupi evaluated for the averages of 
(m^,{m^^}), and similar derivatives for F^'^p^^ and G^'"''\ Equations given by Eqs. (C6)- 
(C13) denote the result of the fourth-order moment method. The second-order moment 
method was applied to a single HH neuron by RT [19] [20], whose result is given by Eqs. 
(C6)-(C10) when we set i = j = 1, -u; = 0, /9i = and Z^^^ = Z^^up = Zup,uq = 0. Equations 
(C6)-(C13) lead to DMA equation given by Eqs. (45)-(52j and (B1)-(B6) if we adopt the re- 
lations as given by Eqs. (A5)-(A7). In particular, in the case of = 1, Eqs. (C6)-(C13) are 
identical with Eqs. (45)-(52) and (B1)-(B6) of DMA if we read ml = C^;j, = 7^,a = Pk,\, 
and = X,,, = y,,,. 
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FIGURES 



FIG. 1. Time courses of (a) fi^, (b) cr^ (= ■\/^v,v) and (c) ag (= ■\/pv,v), for Po = 0.1, /?i = 0, 
J = and N = 100, K^^) and Uo in (a) being shown in arbitrary units. 

FIG. 2. Time courses of (a) Ze and (b) Zg for /?o = 0.1, /?i = 0, J = and N = 100, 

FIG. 3. (a) The /?o dependence and (b) the /?i dependence of Stoe (squares) and Stgg (circles) 
for /?! = in (a) and /3o = 0.1 in (b) with J = and N = 100, filled symbols denoting results in 
DMA and open symbols those in simulations. 

FIG. 4. Log-log plots of 6to£ (squares) and 6tog (circles) against N for (a) Pi = and (b) 
Pi = 0.05 with Pq = 0.1 and J = 0, filled symbols denoting results in DMA and open symbols 
those in simulations. 

FIG. 5. The J dependence of Sto£ (squares) and Stag (circles) for (a) Pi = and (b) Pi = 0.05 
with Pq = 0.1 and N = 100, filled symbols denoting results in DMA and open symbols those in 
simulations. 

FIG. 6. The time course of synchronization ratio 5* for (a) Pq = 0.1, Pi = and J = 100, (b) 
/3o = 0.1, /3i = and J = 200, and (c) Pq = 0.1, Pi = 0.05 and J = 100 with N = 100, solid curve 
denoting results of DMA and dashed curve those of simulations. 

FIG. 7. The dependence of the maximum of the synchronization ratio on (a) J, (b) N, (c) 

Po and (d) Pi: filled and open circles denote Sj^ax 

of DMA and simulations, respectively: filled 

squares express S'^g^^ of DMA (see text) . 

FIG. 8. Time courses of (a) fi^ and (b) 7^,,^ for Poisson spike inputs with the average ISI of 
25 ms for Pq = 0.1, Pi = 0, J = and = 0, solid and dashed curves in (a) denoting results of 
DMA and simulations, respectively. K^^") and Uo in (a) is plotted in arbitrary units. The result of 
simulations in (b) is shifted upwards by 30. 

FIG. 9. Time courses of (a) /i^ and (b) (= \/^v,v) with Pq = 0.1, Pi = and A?^ = 1 for 
constant current input of Ij = 20, solid, dotted and dashed curves denoting results of DMA, DMA2 
(the second-order DMA) and simulations, respectively. A constant input current is shown at the 
bottom of (a) (see text). 

FIG. 10. Time courses of (a) and (b) (= ^/'yv,v) with Pq = 0.2, Pi = and A" = 1 for 
a periodic spike train input with ISI of 25 ms, solid, dotted and dashed curves denoting results 
of DMA, DMA2 (the second-order DMA) and simulations, respectively. A periodic input spike is 
shown at the bottom of (a) (see text). 



21 



This figure "Figl-5.GIF" is available in "GIF" format from: 



http://arXiv.org/ps/cond-mat/0302502v2 



This figure "Fig6-8.GIF" is available in "GIF" format from: 



http://arXiv.org/ps/cond-mat/0302502v2 



This figure "Fig9-10.GIF" is available in "GIF" format from: 



http://arXiv.org/ps/cond-mat/0302502v2 



